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SECTION  1 

INTRODUCTION 


This  program  is  directed  to  a  theoretical  treatment  of 
detonation  behavior  by  numerical  integration  of  the  basic  time  de¬ 
pendent  equations  which  define  the  detonation  process.  Ultimate  ob¬ 
jective  is  to  perform  such  integration  in  two  space  dimensions  so  that 
the  effect  of  geometry  upon  behavior  (for  example  the  phenomena  of 
critical  size)  can  be  predicted  in  the  case  of  large  charge  sizes. 

The  detailed  form  of  the  basic  equations,  and  numerical  procedures 
used  in  their  integration  have  been  described  in  previous  reports.1 


During  the  present  report  period,  one  minor  refinement  has 
been  made  in  the  mathematical  procedures.  This  was  primarily  in  the 
interest  of  completeness,  rather  than  any  anticipated  significant 
effect  on  the  results  of  the  integration.  An  important  development 
has  been  the  discovery  that  the  use  of  the  von  Neuman  "q"  in  inte¬ 
grating  through  a  reactive  shock  can,  under  some  circumstances, 
falsify  the  entropy  increase  and,  as  a  consequence,  seriously  affect 
the  computed  final  steady  state  propagation  velocity.  This  is  now 
regarded  as  the  most  likely  explanation  for  the  high  computed  steady 
state  velocities  discussed  in  the  last  quarterly^  and  previous  re¬ 
ports^.  Some  changes  in  procedure  will  be  required  as  a  consequence 
of  this  discovery. 


The  two  dimensional  program  is  now  in  the  final  checkout  and 
debugging  stage.  It  is  being  actively  pressed  and  should  be  completed 
soon.  This  will  be  affected  somewhat  by  the  problem  associated  with 
the  von  Neuman  "q"  described  above. 

The  equation  of  state  work  is  progressing  at  an  accelerated 
pace.  Data  are  presented  in  this  report  on  polyurethane.  Procurement 
of  data  on  inorganic  salts  such  as  ammonium  perchlorate,  and  other 
organic  binders,  and  composite  mixtures  of  the  two  is  presently  underway. 
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The  differencing  scheme  for  evaluation  of  is  not  of  high  accuracy. 
However,  it  is  regarded  as  suitable  for  present  purposes. 


As  expected,  the  effect  of  this  correction  on  wave  behavior  is 
found  to  be  negligible  for  the  cases  presently  under  treatment.  This  is 
shown  by  the  wave  trajectory  of  Figure  1,  curve  "a"  which  was  computed 
using  the  corrected  program  and  the  same  input  data  as  for  the  trajectory 
shown  on  Figure  5  of  Reference  1.  The  two  curves  are  essentially  identical. 


2.2  STEADY  STATE  PROPAGATION  VELOCITY 


Recent  one  dimensional  calculations  on  TNT  have  encountered 
difficulty  in  that  the  ideal  steady  state  detonation  velocities  obtained 
with  computed  waves  are  about  12  mm./ju  sec,  60  percent  higher  than  the 
values  predicted  from  steady  state  theory.  The  explanation  for  this  dis 
crepancy  has  been  subjected  to  considerable  scrutiny  during  this  report 
period,  since  the  correct  prediction  of  steady  state  values  is  an  im¬ 
portant  test  of  the  validity  of  the  computational  techniques.  As  dis¬ 
cussed  briefly  in  Reference  2,  one  possible  explanation  for  the  high 
steady  state  velocity  is  that  values  used  for  the  equation  of  state 
parameters  are  incorrect. 


The  sensitivity  of  final  wave  propagation  velocity  to  one  of 
the  equation  state  parameters  was  examined  by  comparison  of  the  tra¬ 
jectories  obtained  at  two  difference  values  of  the  parameter  "b"  which 
defines  y  according  to  the  expression 
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FIGURE  1.  EFFECT  OF  PARAMETER  b  ON  WAVE  TRAJECTORY 


c ^yrdj^tocf^xmi/icm^ 

AERONUTRONIC  DIVISION 


These  were  b  =  1.03  (the  value  used  in  all  previous  calculations  on  this 
program)  and  b  =  0.66.  At  the  specific  volumes  appearing  in  the  problems, 
these  values  of  b  correspond  to  maximum  7' s  of  5.3  and  3.8,  respectively. 

Results  of  the  calculation  are  shown  on  Figure  1,  curves  a  and  b. 
Other  input  data  were  as  indicated  in  the  previous  section  for  curve  "a". 

Use  of  b  =  0.66  leads  to  a  steady  state  velocity  of  7.9  mm,//i  sec,  much 
closer  to  the  value  of  about  7  mm./jj,  sec  deduced  from  the  theoretical  steady 
state  treatment  than  the  11.8  mm.//j  sec  shown  for  b  =  1.03.  It  is  clear 
that  a  slight  additional  adjustment  in  "b"  would  produce  the  correct 
behavior. 


The  value  of  b  =  1.03  was  deduced,  as  discussed  in  Reference  2 
from  some  data  on  the  equation  of  state  of  explosive  reaction  products 
obtained  by  M.  A.  Cook.  This  data,  and  in  particular  its  use  for  the 
equation  of  state  of  the  solid  phase  is  not  beyond  question.  Nevertheless, 
before  attributing  the  high  steady  state  velocities  to  uncertainty  in  this 
data  or  its  application  to  the  present  problem,  it  was  thought  desirable 
to  re-examine  the  mathematical  procedures  for  other  possible  sources  of 
error.  It  has  appeared  that  the  steady  state  velocities  may  depend  in  a 
very  significant  way  upon  the  details  of  the  numerical  integration  in  the 
vicinity  of  the  shock  front,  and  that  methods  now  used  may  be  inadequate 
in  this  region. 

The  integration  presently  makes  use  of  a  procedure  devised  by 
von  Neuman  and  Richtmeyer^  for  integrating  numerically  through  the  math¬ 
ematical  discontinuity  which  exists  at  a  shock.  Their  procedure  was  to 
add  a  quantity  "q"  to  the  equations  of  motion.  The  "q"  is  generally 
likened  to  an  artificial  viscosity  and  is  so  constituted  that  it  converts 
the  shock  discontinuity  into  a  continuous  function  which  changes  rapidly 
over  a  few  space  zones  from  the  situation  in  front  of  the  shock  to  the 
situation  behind  the  shock.  It  is  designed  to  produce  the  correct  shock 
energy,  velocity,  and  pressure  relations  before  and  after  the  regions  in 
which  it  is  active,  but  falsifies  these  variables  within  this  region.  It 
was  shown  also  to  provide  for  the  correct  work  loss  due  to  entropy  increase 
across  the  shock  region. 

The  so-called  von  Neuman  "q"  was  originally  used  with  unreactive 
shocks  and  provides  for  the  correct  treatment  of  such  shocks  by  numerical 
methods.  In  the  case  of  reactive  shocks  such  as  occur  in  detonation,  it 
can  again  be  shown  that  the  use  of  the  "q"  in  the  conservation  equations 
does  not  affect  the  energy,  pressure  and  velocity  quantities  on  the  two 
sides  of  the  front.  However,  in  this  case,  definition  of  a  propagation 
velocity  requires  an  additional  relation  generally  known  as  the  Chapman- 
Jouguet  condition.  It  now  appears  that  under  some  circumstances  use  of 
the  von  Neuman  "q"  can  lead  to  an  incorrect  Chapman -Jouguet  condition. 

This  point  can  be  understood  by  means  of  the  following  arguments. 
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Figure  2a  shows  an  actual  pressure  profile  through  a  detona¬ 
tion  wave  as  envisioned  in  terms  of  the  hydrodynamic  theory.  The  wave 
is  initiated  as  a  shock  in  which  the  pressure  rises  discontinuously  from 
Pq  to  p„.  Reactions  start  at  pressure  P2  and  continue  until  completion 
at  the  Chapman-Jouguet  point,  p,. 

Figure  2b  shows  the  same  wave  computed  using  equations  which 
include  the  von  Neuman  "q".  Here  the  wave  initiates  with  a  gradually 
rising  continuous  pressure  front  which  has  a  rounded  off  peak  pres¬ 
sure,  P2 .  In  between  there  is  a  region  of  steeply  rising  pressure. 

The  "q"  used  in  the  calculation  is  computed  continuously  through  the 
rising  pressure  region  as  proportional  to  the  square  of  the  rate  of  com¬ 
pression,  and  provides  for  the  same  total  energy  deposition  at  the  peak 
of  the  wave,  p2 5  as  in  a  shock  process.  In  both  cases  the  propagation 
velocity  is  given  by 


With  the  situation  represented  as  in  Figure  2b,  it  is  evident 
that  if  chemical  reaction  in  the  charge  is  readily  initiated,  such 
initiation  may  take  place  during  the  initial  part  of  the  pressure  rise 
and  reaction  may  be  largely  complete  by  the  time  the  material  has  passed 
through  the  front  and  reached  the  point,  P2.  Thus,  in  contrast  to  the 
situation  represented  in  Figure  2a,  a  major  part  of  the  chemical  reaction 
and  the  corresponding  heat  release  occurs,  in  effect,  as  part  of  the  wave 
shock  process.  This  has  been  observed  to  occur  in  typical  examples  of 
computed  waves . 


The  consequences  of  this  behavior  in  terms  of  wave  velocity 
can  be  appreciated  by  consideration  of  a  typical  Hugoniot  diagram  as 
shown  on  Figure  3.  According  to  the  conventional  interpretation^’^, 
Figure  2a  is  represented  on  such  a  diagram  as  a  process  in  which  an  in¬ 
crement  of  material  starts  at  0  as  undisturbed  material,  jumps  to  2  on 
passing  through  the  lead  shock  and  then  moves  down  the  line  21  while 
undergoing  reaction,  to  1,  the  Chapman-Jouguet  point.  At  point  2, 
energy  deposited  by  the  shock  in  unreacted  material  is  given  by  area  042. 
In  the  subsequent  expansion  to  point  1,  work  represented  by  area  4216  is 
performed  and  the  reaction  energy,  Q,  is  released.  Thus,  at  point  1, 
energy  remaining  in  the  system  is  given  by  area  061  +  Q.  In  the  repre¬ 
sentation  of  the  process  shown  in  Figure  2b  on  a  Hugoniot  diagram,  it 
is  convenient  to  consider  the  extreme  case  in  which  reaction  has  gone  to 
completion  by  the  time  point  2  is  reached.  In  this  circumstance,  points 
1  and  2  coincide  and  on  Figure  3  the  system  goes  directly  from  point  0  to 
point  1.  Since  this  is  treated  as  a  shock  transition,  energy  equivalent 
to  area  061  is  deposited  in  the  process.  The  system  energy  at  point  1 
therefore  is  equivalent  to  area  061  +  Q,  exactly  as  before. 


-7- 


(b) 

FIGURE  2.  WAVE  PROFILES  WITH  AND  WITHOUT  VON  NEUMAN  0 
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FIGURE  3.  HUGONIOT  DIAGRAM  FOR  DETONATION  PROCESS 
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The  first  process  can  be  considered  as  an  isentropic  com¬ 
pression  along  the  adiabat  053,  followed  by  heating  (via  shock  dis¬ 
sipation)  along  32  with  a  corresponding  increase  in  entropy  given  by 

r 

I  TdS  =  area  032,  thence  isentropic  expansion  along  the  adiabat  27, 

J3 

followed  by  heating  due  to  reaction  along  71,  the  latter  involving  a 
further  entropy  increase  given  by  I  TdS  =  Q  (the  heat  of  reaction). 

h 

The  second  process  is  made  up  of  compression  along  the  adiabat  05, 
and  heating  along  51  (due  to  chemical  reaction  and  shock  dissipation) 


with  entropy  increase  corresponding  to  TdS  =  area  051.  Considering 

'  5 

that  temperatures  are  highest  at  the  Chapman-Jouguet  point,  it  is 
apparent  from  these  relations  that  the  entropy  increase  in  the  second 
process  is  smaller  than  the  first.  Since  the  final  energy  content  is 
the  same  in  both  cases,  it  must  be  concluded  that  the  second  process 
results  in  a  higher  pressure  at  the  Chapman-Jouguet  point  than  does  the 
first.  The  second  process  cannot,  therefore,  terminate  at  point  1,  but 
must  instead  arrive  at  some  point  at  a  higher  pressure,  such  as  8.  From 
equation  (1) ,  the  relative  propagation  velocities  in  the  two  cases  will 
be  given  by  ratio  of  the  square  root  of  the  slopes  of  the  lines  01  and  08. 
It  thus  becomes  evident  that  the  situation  represented  by  Figure  2b  will 
correspond  to  a  higher  velocity  than  the  value  deduced  from  the  steady 
state  treatment. 


The  foregoing  considerations  further  suggest  that  propagation 
velocities  as  deduced  from  the  steady  state  treatment  will  be  obtained 
if  reaction  is  not  permitted  to  occur  until  the  system  has  passed  com¬ 
pletely  through  the  shock  zone,  i.e.,  until  point  p£,  Figure  2b,  or 
point  2,  Figure  3,  is  reached.  In  effect,  this  forces  the  system  to  go 
through  the  sequence,  shock  compression  from  point  0  to  1  (Figure  3)  and 
expansion  with  reaction  along  line  21,  as  is  implied  in  the  steady  state 
treatment.  This  has  been  subjected  to  verification  by  computing  a  wave 
trajectory  in  which  the  IBM  709  was  arbitrarily  programmed  to  provide 
such  a  reaction  delay.  Results  are  shown  on  Figure  4  as  curve  1.  For 
comparison  purposes,  curve  2  is  included^which  is  the  trajectory  shown 
on  Figure  5  of  the  last  quarterly  report  ,  obtained  with  no  restrictions 


on  reaction.  Other  input  data  was  identical  for  the  two  curves  and  is  as 
previously  reported^. 


The  difference  between  the  two  curves  is  pronounced,  and  in¬ 
dicates  that  the  effect  discussed  is  very  important.  Because  of  the  long 
reaction  zone  existing  in  the  wave,  curve  1  never  accelerated  to  steady 
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state.  It  is  therefore  not  ascertained  whether  the  programmed  reaction 
delay  leads  to  the  correct  steady  state  velocity.  However,  it  is  clear 
that  the  details  of  the  mathematical  treatment  in  the  vicinity  of  the 
wave  shock  zone  are  more  important  than  has  heretofore  been  thought  to 
be  the  case.  A  number  of  procedures  for  obtaining  a  more  correct  treat¬ 
ment  are  possible  and  these  are  presently  being  considered. 

It  should  bet  noted  that  this  problem  does  not  affect  the 
validity  of  the  calculations  on  minimum  pressure  for  initiation  of  ex¬ 
plosive  materials,  as  discussed  in  previous  reports^’  .  This  is  because 
in  those  cases,  the  critical  point  was  that  at  which  reaction  ceased  in 
the  wave.  At  this  point,  the  wave  becomes  an  unreactive  shock,  and  the 
present  procedures  are  correct  for  an  unreactive  shock. 
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SECTION  3 


TWO  DIMENSIONAL  PROGRAM 


Mathematics  and  coding  are  essentially  complete  and  de¬ 
bugging  is  now  in  process  on  this  program.  However,  ^.t  is  subject  to 
the  considerations  discussed  in  the  previous  section  and  may  have  to 
be  delayed  until  improved  procedures  for  treating  the  shock  transition 
are  developed.  Any  such  delay  should  be  of  short  duration.  Every 
effort  is  being  made  to  get  this  program  running  as  soon  as  possible. 
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SECTION  4 

EQUATION  OF  STATE  STUDIES 


The  equation  of  state  studies  have  involved  some  work  on 
the  refinement  of  techniques  and  calibration  of  equipment.  However, 
a  major  part  of  the  effort  during  this  report  period  was  directed  to 
studies  on  the  equation  of  state  of  polyurethane. 

Methods  were  as  described  in  the  last  quarterly  report. 

Data  is  obtained  in  the  form  of  a  plot  of  pressure  vs.  piston  movement 
over  the  up  and  down  legs  of  a  pressurization  cycle.  These  legs  diverge 
from  the  ambient  pressure  point  due  to  a  combination  of  piston  sliding 
friction,  and  sample  internal  friction.  This  is  the  same  in  both  di¬ 
rections^  and  therefore  symmetrical  around  a  center  curve  except  for  a 
distorted  region  at  each  end  due  to  the  reversal  of  motion. 

The  magnitude  of  the  frictional  forces  is  indicated  by  the 
horizontal  separation  of  the  up  and  down  legs  in  the  symmetrical  region. 
Press  sliding  friction  is  a  function  of  pressure  but  independent  of 
specimen  material  so  that  variations  in  separation  for  different  mater¬ 
ials  or  conditions  are  normally  due  to  differences  in  sample  internal 
friction,  i.e.,  viscosity.  Information  on  this  property  is  of  signif¬ 
icant  scientific  interest,  particularly  in  the  present  instance  since  it 
may  afford  some  insight  into  possible  differences  between  the  equation 
of  state  observed  under  the  present  static  methods,  and  that  which  applies 
to  conditions  of  shock  loading  as  encountered  in  a  detonation  process. 

No  attempt  at  a  detailed  analysis  on  this  basis  has  been  made.  The 
question  will  be  examined  further  at  a  later  date. 

Data*on  two  typical  press  cycles  to  26  and  40  kilobars  and  at 
ambient  temperature  are  shown  in  Figure  5  as  plots  of 1/1  vs.  pressure 
(1  ==  sample  length).  The  values  of^l  include  a  correction  for  press 


*Data  on  the  polyurethane  specimen  are  shown  in  Table  I. 
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TABLE  I 

PREPARATION  OF  POLYURETHANE  SPECIMEN 


Polymer 


Polypropylene  glycol  2025 

83.30% 

Trimethylol  propane 

3.19% 

Ferric  acetylacetonate 

0.125% 

Toluene  -  2,4  -  diisocyanate 

13.85% 

Cure 

24  hours  @  160  F 

Sample 


Density 

1.068  g/cm3 

Sample  wt. 

5.0532  g 

Sample  initial  volume 

4.729  cm3 
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FIGURE  5.  PRESS  CYCLES  FOR  POLYURETHANE  TEMP.  -  23.3°C  0.5  INCH  CYLINDER 
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distortion  obtained  by  running  a  similar  cycle  with  no  specimen  in  the 
pressure  chamber.  The  midpoint  between  the  up  and  down  legs  is  indi¬ 
cated  as  curves  231  and  23H  which  are  extrapolated  (dashed  portion)  be¬ 
yond  the  region  of  symmetry.  The  validity  of  the  symmetry  assumption 
is  supported  by  the  close  agreement  of  these  centered  curves. 

The  press  cycle  and  centered  curve  for  polyurethane  at  100°C 
is  shown  in  Figure  6.  Evident  here  is  an  increase  in  compressibility  and 
a  decrease  in  internal  friction,  the  latter  shown  by  the  reduced  separa¬ 
tion  between  the  up  and  down  legs.  Increase  in  compressibility  with 
temperature  is  normal  behavior  for  many  organic  substances. 

Final  p-v  curves  for  the  ambient  and  100°  temperatures  are 
shown  in  Figure  7.  A  slight  increase  in  volume  as  well  as  the  increase 
in  compressibility  with  temperature  is  noted  at  ambient  pressure.  Be¬ 
cause  of  this,  the  curves  for  the  two  temperatures  cross  at  about  2.5  kbar. 
In  the  region  of  25  to  40  kbar  the  compressibility  decreases  with  temper¬ 
ature  as  would  be  expected. 

Attempts  at  interpretation  of  this  work  in  terms  of  a  general 
equation  of  state  will  be  deferred  until  more  complete  data  becomes 
available. 
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FIGURE  6.  PRESS  CYCLE  FOR  POLYURETHANE.  MAXIMUM  PRESSURE  - 

40  KBAR,  TEMP.  RANGE  99.5  -  101°C  0.5  INCH  CYLINDER 
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FIGURE  7.  PRESSURE  -  VOLUME  CURVES  FOR  POLYURETHANE 
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